Sweeping from the superfluid to Mott phase in the Bose-Hubbard model 
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We study the sweep through the quantum phase transition from the superfluid to the Mott state 
for the Bose-Hubbard model with a time-dependent tunneling rate J(t). In the experimentally 
relevant case of exponential decay, J(t) oc e -7 ', an adapted mean-field expansion for large fillings n 
yields a scaling solution for the fluctuations. This enables us to analytically calculate the evolution of 
the number and phase variations (on-site) and correlations (off-site) for slow (7 <C fi), intermediate, 
and fast (non-adiabatic 7 3> fi) sweeps, where /i is the chemical potential. Finally, we derive 
the dynamical decay of the off-diagonal long-range order as well as the temporal shrinkage of the 
superfluid fraction in a persistent ring-current setup. 
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For many systems, the equilibrium properties (the 
thermal ensemble at finite temperatures or the ground 
state at zero temperature) are quite well understood. 
However, strictly speaking, the equilibrium state applies 
to purely static situations only. For dynamical systems, 
the adiabatic theorem states that the actual quantum 
state remains close to the ground state (at zero temper- 
ature) if the external time-dependence is slow enough, 
i.e., much slower than the internal frequencies of the sys- 
tem determined by its energy gaps. If this adiabatic- 
ity condition fails, however, non-equilibrium effects be- 
come important, leading to many intriguing phenomena. 
For example, the usual split into the ground state plus 
small quantum fluctuations around it is no longer unique, 
which usually leads to effects such as the amplification of 
quantum fluctuations and the creation of quasi-particles. 

A prototypical example for such a situation is a second- 
order quantum phase transition Q, cf. Fig.Q] At the 
critical point, the energy gap vanishes and the response 
time diverges. Consequently, sweeping through the phase 
transition by means of a time-dependent external param- 
eter 1/J(t) with a finite velocity dJ/dt < inevitably 
violates adiabaticity close enough to the critical point 
and so generates non-equilibrium effects. Fostered by 
the tremendous progress in the experimental capabilities, 
there has been increasing interest in quantum phenomena 
at low temperatures in general and quantum criticality in 
particular. In view of the creation of entanglement 
quantum phase transitions are also relevant for quan- 
tum information (the continuum limit of an adiabatic 
quantum algorithm represents just a sweep through a 
quantum phase transition). Finally, dynamical quantum 
phase transitions bear strong similarities to cosmological 
phenomena: If the stable phase permits topological de- 
fects (symmetry-breaking transition) , a quantum version 
of the Kibble-Zurek mechanism may create those defects 
0,0, 00| In case the initial spectrum contains gapless 
modes, their quantum fluctuations will generally be am- 
plified in close analogy to cosmic inflation 0, 0] . Both 



phenomena can be understood in terms of the emergence 
of an effective cosmic horizon which corresponds to the 
loss of causal connection and the breakdown of adiabatic- 
ity near the critical point. 
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FIG. 1: [Color online] Sketch (not to scale) of level structure 
near the critical point, i.e., plot of the energy of the ground 
state (green solid line) and the energy gap (thin blue line) as a 
function of the external parameter 1/J(t). Since the response 
time diverges at the critical point J = J c , adiabaticity breaks 
down near the transition for any finite sweep velocity dJ/dt. 
Hence the actual state j*(t)) (red dashed line) slightly devi- 
ates from the exact ground state already before the critical 
point J = J c is reached and finally lies energetically some- 
where between the coherent state (dotted green line) and the 
Mott state. 

Besides the Ising model (which can be solved analyti- 
cally 0,0]), another archetypical example is the quantum 
phase transition between the superfluid and the Mott 
insulator state in the Bose-Hubbard model (where no 
closed solution has been found). This model provides 
a simplified description for the problem of interacting 
bosons hopping on a lattice ; its implementation with 
atoms/molecules in optical lattices [12], and the sub- 
sequent experimental realization of the Mott-Hubbard 
transition [lj| has recently caused a large amount of re- 
search activity ^3] • The superfluid phase possesses a rich 
spectrum ranging from topological defects (e.g., kinks 
in one dimension) to gapless quasi-particle excitations 
(phonons) but the Mott state does not. Hence the (quan- 
tum) Kibble-Zurek mechanism requires the dynamical 
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transition from the Mott to the superfluid state (quan- 
tum melting) whereas the amplification of the quantum 
fluctuations of the phonon modes occurs in the opposite 
direction (quantum freezing). The (mostly numerical) 
studies thus far 0, have been devoted to the first 
case (quantum melting). Here, we study the dynamical 
behavior of the superfmid to Mott quantum phase tran- 
sition starting, by contrast to previous investigations, 
from the superfluid side of the transition. We thus in- 
vestigate the quantum freezing of number fluctuations 
[la,li3 an d predict their dependence on the sweep rate 7. 
The thereby created number-squeezed state and its num- 
ber fluctuations can be studied experimentally, for ex- 
ample by the interference of condensates released from 
different optical lattice wells As further observables 
representing the actual quantum state \^(t)}, to bring 
out the difference to either the superfluid or Mott states 
(cf. Fig.^J, we derive the temporal evolution of the off- 
diagonal long-range order and the superfluid fraction. 

The Bose-Hubbard model with local contact interac- 
tions is, on a rather general type of lattice, described by 
the Hamiltonian 



H = J(t) £ M a0 al a? + |$>t) 2 a 2 



(1) 
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Here and a a denote the bosonic creation and an- 
nihilation operators at the lattice site a. The energy 
scale is set by the on-site repulsion U > and the time- 
dependent tunneling rate J{t) > 0. The structure of the 
lattice is encoded in the matrix M a p, which is supposed 
to obey the same (translational) symmetry group, e.g., 
for a one-dimensional chain with nearest-neighbor hop- 
ping, we have M al3 = 5 a ,p - \{8 a ,f}+i + S a ,p-i). 

The initial state is the homogeneous superfluid phase, 
J S> J c — 0(U/n), see [lU, with a large integer filling 
n = (n a ) ^> 1 (where h a — aj,a a ), which facilitates a 
particle- number conserving mean-field expansion (511 Ell 



(Vo+Xa+Ca)^/^, 



(2) 



with A = a s (a^a s )- 1 / 2 7V 1 / 2 and iU = N = £ Q n Q , 
where as = J2a®a- Here the main idea is to expand 
the original operator a a into powers of n ^ 1, consider- 
ing the (formal) limit n | 00 with the chemical potential 
[i = Un remaining finite, i.e., U — 0(l/n). The lead- 
ing term is the order parameter i/'o = ^(y/n) and the 
quantum corrections Xa — O(n ) correspond to quasi- 
particle excitations. For the validity of this expansion, 
the remaining higher-order corrections £ Q must be small 
Cq <C O(n ). Inserting the above expansion Q into the 
equation of motion ihdta a = J(t) M a pdp + Un a a a 
obtained from Q and sorting into powers of n yields the 
analogue of the Bogoliubov-de Gennes equations for the 
excitations at site a (h = 1 throughout) 

idtXa = J(t) M ^XP + 2C/|^ 2 |xa + U^lxl , (3) 

P 



and for the remaining higher-order corrections 

id t ( a = J(t)Y,M a pCp + 2U\^\( a + U^ 2 (i+ (4) 



+w%xix a + ur xl + uxixl + 0(U( a ) • 

Deep in the superfluid phase (our initial state), the 
higher-order corrections £ Q are small and the mean-field 
expansion l|2| works very well. If we approach the Mott 
phase, however, these corrections start to grow accord- 
ing to Eq. Q and, at some point, the mean-field ex- 
pansion (2) breaks down. The characteristic time-scale 
of this breakdown can be estimated from the nonlinear 
source terms in Eq. <@J which are suppressed to 0(l/y/n) 
in view of U = 0(l/n), = 0(^/n) and x a = O(n ). 
Hence (starting in the superfluid phase), the higher-order 
corrections remain small as long as Uty/n-^L 1, i.e., for 
evolution times of order t = 0{\fn). Thus we can extra- 
polate the mean-field expansion even into the Mott phase 
for some time t = 0(*/n), and follow the evolution of the 
instabilities which develop because the superfluid state is 
no longer the ground state of the system. 

The polar decomposition a a — exp{i(f) a }\/n^ yields 
the linearized number fluctuations 6n a via h a = n + Sn a 
and accordingly the conjugate phase fluctuations 5(j) a , 
where x Q = ipo[Sn a /(2n) + i5<j) a ] + 0(l/y/n). In terms of 
these fluctuations Sh a , Eq. (EJ can be diagonalized by a 
normal-mode expansion into the eigenvectors of the ma- 
trix M a p labeled by the generalized momenta K 



5n K 



0, (5) 



where X K are the corresponding eigenvalues of M a p. 

For small A K , the above evolution equation is analogous 
to the modes of a quantum field within an expanding 
universe with the wavenumber k oc \/Ak' Pursuing the 
similarity a bit further, we get the analogue of a cosmic 
horizon with the horizon size 



A h (t) = / dt' y/J(t')Un. 



(6) 



if J(t) changes fast enough 0, 0|. Note that the up- 
per integral limit should not be taken literally, since it 
must still be within the region of validity of the mean- 
field expansion. In case an effective horizon occurs, its 
size constantly decreases. Hence all modes (with small 
\ K ) will qualitatively follow the same evolution - but at 
different times: Initially, the wavelength of the modes is 
well inside the horizon, AkA^ ^> 1, and the modes os- 
cillate almost freely. At some point of time, however, 
the constantly shrinking horizon closes in, A re A 2 = 0(1), 
and thus the causal connection across a wavelength is 
lost. After that, the modes cannot oscillate anymore and 
freeze. This process leads to the amplification of the ini- 
tial quantum fluctuations via squeezing in analogy to the 
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inflationary epoch in the early universe (SlLLJJ]. We note 
that these considerations show a major weakness of the 
adiabatic-impulse approximation used in Q , since modes 
with different X K become non-adiabatic at different times. 

Since the tunneling rate J depends exponentially on 
the amplitude of the laser which generates the opti- 
cal lattice (52], we shall consider an exponential time- 
dependence J(t) = Joexp{— 7^} in the following, which 
is implying the emergence of a quasiparticle horizon ac- 
cording to Eq. Jsj. Furthermore, Eq. (§J reveals that X K 
can be absorbed by a suitable redefinition of the time 
coordinate in this — AX R J(t)/j, leading to 

( d 2 |\ 2 Un\\ _ A 
TT + 1 )Sn K = 0. 7 

The universal scaling solution resulting from the above 
equation confirms that modes with different A K display 
the same behavior, but at different times. The only re- 
maining dimensionless parameter is the ratio v = Un/j 
which is a measure of the sweep velocity: v 3> 1 implies 
a slow and v -C 1 a fast (non-adiabatic) sweep. The dif- 
ferential equation <|2jl has a universal scaling solution in 
terms of Whittaker functions ji^] 

Sh K = Vn~e~™ /2 W^ 1/2 (2iT K ) b K + h.c. (8) 

The quasi-particle operator b K annihilates the adiabatic 
ground state b K |in) =0 at early times t k J, — oo, where 
the modes oscillate like e ±iTK (which can also be ob- 
tained from the asymptotic behavior of the Whittaker 
functions). According to the arguments presented after 
Eq. JIJ , the mean-field expansion remains valid for inter- 
mediate times with Uty/n <C 1 but jt^$>l, and hence can 
be extrapolated to the "late-time" regime t k f where 
the functions W il/i i/ 2 (2zT K ) approach a constant value, 
i.e., the number fluctuations freeze (J J. 0). Due to the 
perfect scaling solution, the frozen value is independent 
of k, but the decaying corrections do depend on X K : 

(Shi) = H (Sh K ) 2 |in) = n 1 " 6 2 ™ + OitX^e-' 1 ) . (9) 

Since the leading term is independent of k, it just yields a 
local (oc § at p) contribution after the mode sum (k — > a) 
and thus leads to frozen on-site number variations 

A 2 (n Q ) = (hi) (h a ) 2 = (8hl) = n l —f^ . (10) 

For a rapid sweep v <C 1, this variation A 2 (n a ) ap- 
proaches a constant value n, which is characteristic of a 
superfluid phase with Poissonian number statistics 
For a slow sweep v ^> 1, the variation A 2 (n a ) becomes 
small, oc n/v, and so approaches the behavior deep into 
the Mott phase (with J | 0), where A 2 (n a ) = 0, cf. [H]. 
For J | 0, the energy (H) is basically determined by the 
variation A 2 (n a ) and hence this quantity describes the 



energetic location of the final state in comparison to the 
coherent state and the Mott state, see Fig.Q] 

The sub-leading corrections in (8h 2 K ) depend on k and 
hence determine the off-site (a ^ j3) number correlations 
which decay exponentially for ^> 1 

(h a hp) - (n a )(fip) = (8h a 5hp) = O^te^ 1 ) , (11) 

where we have omitted finite-size effects (which scale with 
the inverse number of lattice sites) , since strictly speaking 
the mode sum does not include the zero-mode k = 0. 

The conjugate phase fluctuations can be derived in an 
analogous manner and are determined by the derivatives 
of the Whittaker functions dW^ i^/dr^. As one would 
expect they do not freeze - but increase: 

(S4>1) - v 1 7 2 t 2 +0( 7 tlnA K ) . (12) 

Again the leading (first) term is independent of n and 
thus yields the on-site fluctuations (S(p1) only (which 
generate the quadratically growing quantum depletion 
(XaXa) and anomalous term (Xa))- The second term 
is the leading K-dependent contribution and determines 
the off-site phase correlations (8<j) a 8^)p). 

The ascent of the phase fluctuations can be understood 
as a consequence of the emergence of an effective hori- 
zon which entails the loss of causal connection between 
different sites and thus the decay of the phase coher- 
ence across the lattice. As one interesting observable, 
let us discuss the evolution of the off-diagonal long-range 
order between sites a and j3 defined by the correlator 
(oJ,(t)d/3(<)). Using the mean-field results valid at inter- 
mediate times with -ft ^> 1 but Uty/ri < 1 as the initial 
conditions, we may derive the ensuing stages of the quan- 
tum evolution, where the tunneling rate J(t 3> I/7) <ggC 1 
is exponentially small and can be neglected. In this 
limit, the evolution of the operators can be approxi- 
mated by da a /dt = —iUh a d a which possesses the sim- 
ple solution d a (t) = exp{— iUh° a t}d a a . Consequently, we 
obtain (a£(i)a0(t)) = n(exp{iU(h a — hp)t}) + 0(y/n). 
The frozen first-order number fluctuations 8h a are in a 
squeezed state which can (for n ^> 1) be approximated 
by a (continuous) Gaussian distribution. For a Gaussian 
variable X with (X) = 0, the exponential average yields 
(cxpjzA}) = cxp{ — (A 2 )/2} and hence we get 

(al(t)ap(t)) pb nexp{-U 2 t 2 A 2 (n a )} . (13) 

Note that this expression is only valid for time scales 
I/7 < t < 1/U; for Ut € 27rN all the exponential factors 
equal unity again and we have a revival of the phase co- 
herence and thus long-range order similar to a spin echo. 
(Note that we omit the coupling to the environment in 
our derivation.) Furthermore, the above Gaussian decay 
is independent of the distance between the sites a and /3. 
Since the Fourier transform of determines 
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the structure factor S(k), the decay of the long-range or- 
der (1 1 -3 1 directly corresponds to the temporal decrease of 
the peak in S(k) at k — 0. This dependence on time and 
sweep rate can be observed with time-of-flight measure- 
ments [26] (which basically map out |5(fe)|) by varying 7 
and the time delay between the phase transition and the 
release of the condensate. 

Let us calculate another simple experimental signa- 
ture, which involves just neighboring sites, namely the 
superfiuid fraction n a f/n. To this end, we specify our lat- 
tice and choose a one-dimensional rin g-ge ometry a — > I 
with the circumference coordinate I [22]]. In order to 
generate a persistent current, we impose a small phase 
gradient externally by trapping a flux quantum on the 
ring: = y/nexp{—i[it + 2iri£/L}, where L ^> 1 de- 
notes the number of sites. Naturally, for a decreasing 
tunneling rate J(t), the induced current diminishes also. 
In addition, the number of particles contributing to this 
flux goes down - which will be used as a measure for the 
superfiuid fraction. The flux can be determined by the 
usual Noether current corresponding to the U(l) invari- 
ance and is related to the nearest neighbor correlation 
function (aja^+i). Hence, the shrinkage of the superfiuid 
fraction is given by the same j3| expression as in Eq. l|T3jl 

T"- p {- oV - i ^}- <14 » 

In the case of a very rapid (impulse) sweep, i/<l, the 
superfiuid fraction decays with exp{— nU 2 t 2 }, i.e., inde- 
pendent of 7. Conversely, for a slow sweep v ^> 1, the 
decay takes much longer: exp{— nU 2 t 2 / (2irv)} . 

In conclusion, an adapted mean-field expansion en- 
ables us to analytically calculate the freezing of number 
fluctuations and the growth of phase fluctuations during 
an (exponential) sweep through the superfiuid — > Mott 
quantum phase transition. As further experimental sig- 
natures, we predict the temporal decay of the superfiuid 
fraction n s f/n and of the central k = peak in the struc- 
ture factor S(k). 
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